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FOREWORD 


This report describes work performed by Systems and Applied 
Sciences Corporation (SAS) under NASA Contract No. NAS5-24255. The work 
involves development of an improved two-dimensional model of chemical 
and physical processes in the stratosphere. 
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SECTION 1 - INTRODUCTION 


A more comprehensive model for stratospheric chemistry and transport 
theory is being developed for the purpose of aiding predictions of changes in the 
stratospheric ozone content as a consequence of natural and anthropogenic 
processes. This new and more advanced stratospheric model is time dependent 
and the dependent variables are zonal means of the relevant meteorological 
quantities which are functions of latitude and height. The model is constructed 
by the best mathematical approach on a large IBM S360 in American National 
Standard FORTRAN. It will be both a scientific tool and an assessment device 
used to evaluate other models. Program blocks are kept under maximum of 
100 statements per routine and each routine has but one function. The report 
sets out the detailed formulation of a numerical model both in physics and 
mathematics. The interactions of dynamics, photochemistry and radiation 
in the stratosphere can be governed by a set of fundamental dynamical equations 
which wiU be described in Section 2. Some physical ideal is obtained from an 
existing FORTRAN computer source deck of Crutzen's two-dimensional model 
of the stratosphere which is subsequently applied in building new and more 
coraprohensive models. Section 3 describes the numerical method used in the 
integration. Some preliminary results and discussions are presented in Section 4. 
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SECTION 2 -• BASIC EQUATIONS 


The model is an integration of an extensive set of fbndamental equations 
which can be used to predict the future state of the atmospheric circulation from 
knowledge of its present state. In an (x, y,^) coordinate system in which x = 
distance eastwards, y = distance northwards, and ^ = log (Pq/p); with pQ = 1000 mb 
and p - pressure, u, v, and w are the three corresponding velocity components 
respectively. In this log-pressure system the eastward component of the 
momentum equation is: 

^ u + u u + V u + w u — tan cp + 2 v sin ip — (1) 

at ax ay a^ a ^ ax 

where a = radius of the earth, f2 = rate of rotation of the earth, <p - latitude, 

Z is geopotential. The thermodynamic equation is: 

-^P + u V w (2) 

at ax ay 

k 

where 0 = T (P(/p) , T = temperature, q is the rate of change of potential 
temperature, K is the ratio of gas constant for dry air R to specific heat of air 
at constant pressure C . The equation for conservation of matter is: 




Since our aim is to deal explicitly with zonal mean quantities only. 
Accordingly let a over bar denote the zonal mean and a dash denote the departure 
therefrom, then equations (1) to (3) yield, after taking an average around a 
latitude circle: 

(T e”^cosQ7) +-^(V r) + -^(W T ) = ** ) - -^(W' rO (4) 

ht Sy Sy 

0 e"^ cos<p) + ) + -^(W 0 ) = q e"^ cos ■|^(V7T?) ~ ^ (W^”) 

( 5 ) 

^ V* + “W = 0 (6) 

oy 

where T = (fta cos ^ + u) a cos <p 

V = V e ^ cos <p 
W = w e ^ cos (p 

Note that the continuity equation (3) takes the simfde form (6) which is the 
advantage of using the perssure coordinates. 

Equations (4) to (6) may become a complete set by adding the y 
component momentum and hydrostatic equations. For our own purpose of not 
introducing the gravity waves we introduce an approximation to the y component 
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momentum equation that the zonal mean wind is in geostrophic equilibrium with 
the zonal mean geopotential field. Combining this with the hydrostatic equation 
leads to the thermal wind relationsliip: 


A. 


f = R exo (-k ^ ) 
20 tan<p 


ay 


( 7 ) 


Equations (4) to (7) can be regarded as a closed set which completely 
determines the relationships among the dependent variables T, 0, V, W. The 
other variables can be written in: terms of the relevant quantities and are held 
as functions of time, latitude and altitude. 
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SECTION 3 - NUMERICAL METHOD 


It is convenient to make the equations dimensionless. The scaling quantities 
is given a subscript s. Equations (4) to (7) become: 

(r e“^slcos<;^>) + ^(T V) + A (8) 



™ (0 V) + ( 0 \V) = q e^^s^cos <;£?+ Q 

d(p hi 




T-~ 


e 

tan^ 


de 

d<p 


( 9 ) 


( 10 ) 


v= - ^ 


w= 

d<P 


ill) 


where H and Q are the divergence of eddy momentum and eddy heat flux, 
respectively. is a stream function which is introduced to rewrite equation (6) 
in the form of equation (11). The method of solution is to eliminate the time 
derivatives from (8) and (9) using (10). Substituting for V and W from (11) 
gives an equation for ijf . 








Q(p e + sin^ <p) 

sin^ <p 


+ 


?. e“*^^^s0^ytan</?] H + 


’s 




^(p . '=’8 


J I^G ^^Scos (p ^fp'^ ^ /tan <p 

OF poo®* 


( 12 ) 
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The ^ and cp used as subscripts denote differentiation with respect to that 
quantity. It may be solved with suitable boundary conditions for ^ which may 
be used in (8) and (9) to get new values of T,0 at the advanced time step. 

Equation (12) is then solved at the new time value and the process repeated. 

A Crank-Nicholson implicit finite difference scheme is used in the 
numerical calculations . The Alternate Direction Implicit method is adopted so 
that the finite difference equation can be solved by simply inverting a tridiagonal 
matrix in each of the two dimensions sequentially. The nonlinear coefficients 
in the equations can be taken into account by using a predictor - corrector 
procedure. A more detailed description is followed. 

Eq. (12) is an elliptic type, second order partial differential equation which 
can be solved by running a marching method for the time-dependent equation out 
to a sufficiently large tim“. Equations (8) and (9) are time-dependent first order 
partial differential equations. A general form for the set of equations can be 
written as: 

®t “ ® ^ + D + E + F (13) 

The ADI scheme for equation (13) is the foll'^wing: 

S s“ + At [a 6| ( V B S^0 + C 80 (S“ ) 

+ D 6^ (S""^^) + E 60 (S^) + f] (14) 
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( 15 ) 


S“+l - s" t At [ A (S“*2 ) + B * C 6 ^ (S“*l) 

+ D 6j (S“*^ ) + E 6^ (S“** ) + f] 

where? ^ 

®jk = “ 2 )/A|^ 

®jk ^ <®j.k+l - 2 %], + )/A0^ 

“ Vl,k )/2A^ 

®jk "<®j,k+l "®j,k-l )/2A<^ 

A simple tridiagonal matrix inversion wttl solve equations (14) and (15) to give 
the solution of equation (13). 
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SECTION 4 - PRELIMINARY RESULTS AND DISCUSSIONS 

Equation (12) is an elliptic 15 ^ 50 , boundary value problem so that the crossing 
termxl/^^ is not a dominant term which can be approximated at the old time step 
as shown in equation (14) and (15). This enables us to use the Crank-Nicholson 
-cheme and the ADI method to solve equation (12) without changing the coordinates 
for getting rid of the crossing term. Equation (12) has been solved without over 
100 total time steps before it reached the steady state solution. 

We have taken the average for each two consecutive values of i// in the 
process of getting fast convergence. The values of xf/ are then used to determine 
V and W which are substituted in equations (8) and (9) for solving r and 9 at the 
advanced time step. This process has been repeated many times. 

The terms in the right hand side of equation (12) are source fimctions which 
drive the rate of heating due to absorption of solar radiation and the rate of cooling 
by thermal radiation. Eddy fluxes of heat and momentum are treated by large- 
scale diffusion coefficieiits. For these fixed driving terms the solutions of the 
set of equations do not converge very well. It may be due primarily to the poor 
initial data available. Therefore, we have to find the steady state solution first 
before seeking any time-dependent solutions. 
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